knitr::opts_chunk$set(echo = FALSE)
if(!require(dplyr)){install.packages("dplyr")}
## Loading required package: dplyr
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
if(!require(readxl)){install.packages("readxl")}
## Loading required package: readxl
## Warning: package 'readxl' was built under R version 4.0.3
if(!require(ggplot2)){install.packages("ggplot2")}
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 4.0.3
if(!require(tidyr)){install.packages("tidyr")}
## Loading required package: tidyr
## Warning: package 'tidyr' was built under R version 4.0.3
if(!require(plotly)){install.packages("plotly")}
## Loading required package: plotly
## Warning: package 'plotly' was built under R version 4.0.3
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
if(!require(sjPlot)){install.packages("sjPlot")}
## Loading required package: sjPlot
## Warning: package 'sjPlot' was built under R version 4.0.3
## Learn more about sjPlot with 'browseVignettes("sjPlot")'.
if(!require(lmtest)){install.packages("lmtest")}
## Loading required package: lmtest
## Warning: package 'lmtest' was built under R version 4.0.3
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 4.0.3
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
if(!require(plm)){install.packages("plm")}
## Loading required package: plm
## Warning: package 'plm' was built under R version 4.0.3
## 
## Attaching package: 'plm'
## The following objects are masked from 'package:dplyr':
## 
##     between, lead
if(!require(broom)){install.packages("broom")}
## Loading required package: broom
## Warning: package 'broom' was built under R version 4.0.3
if(!require(kableExtra)){install.packages("kableExtra")}
## Loading required package: kableExtra
## Warning: package 'kableExtra' was built under R version 4.0.3
## 
## Attaching package: 'kableExtra'
## The following object is masked from 'package:dplyr':
## 
##     group_rows
if(!require(knitr)){install.packages("knitr")}
## Loading required package: knitr
## Warning: package 'knitr' was built under R version 4.0.3
if(!require(gplots)){install.packages("gplots")}
## Loading required package: gplots
## Warning: package 'gplots' was built under R version 4.0.3
## 
## Attaching package: 'gplots'
## The following object is masked from 'package:stats':
## 
##     lowess
if(!require(DiagrammeR)){install.packages("DiagrammeR")}
## Loading required package: DiagrammeR
## Warning: package 'DiagrammeR' was built under R version 4.0.3
if(!require(tableone)){install.packages("tableone")}
## Loading required package: tableone
## Warning: package 'tableone' was built under R version 4.0.3

Contenidos

logo

  • Configuración de datos
  • EXTRA: ggplot2

dplyr y ggplot

  • TIDYVERSE(DPLYR):
    • Manipulación de datos. Similitud a PANDAS (Python)
    • Flujo de trabajo
    • Permite mayor simplicidad en la codificación
    • Concatenación (pipe)
    • Sintaxis entendible
    • Similar a manipulación vía SQL
    • Desventaja: mucho uso de memoria.
  • GGPLOT:

Ejemplo integración

Procesos: bajar base de datos (1), seleccionar columnas (dplyr::select) (2), formatear fecha (desde formato de texto día-mes-año a “Ymd”) (3), filtrar regiones (sólo Coquimbo y Valparaíso) (4), generar figura de serie de datos (Fecha como eje X y Nuevo Caso Confirmado como eje Y), discriminando por Región por colores (5), generamos un área gris del rango intercuartil del resto de regiones, excluyendo a la región metropolitana (al 80% o fill=grey70) con transparencia para que no se superponga (alpha=.4) (6), definimos los colores de las líneas: rojo para Antofagasta, verde para Atacama, azul para Tarapacá y negro para la mediana del resto de regiones (scale_color_manual) (7).

covid19_chile_coq_val <- 
    readr::read_csv("https://raw.githubusercontent.com/ivanMSC/COVID19_Chile/master/covid19_chile.csv") %>% 
    dplyr::select(Fecha, Region, `Nuevo Confirmado`) %>%
    dplyr::mutate(Fecha=as.Date.character(Fecha,tryFormats = "%d-%m-%Y")) %>% 
    dplyr::filter(Region %in% c("Valparaíso","Coquimbo")) %>% 
    ggplot()+
      geom_line(aes(x=as.Date(Fecha), y=`Nuevo Confirmado`, color=Region))+
      scale_x_date()+
      xlab("Fecha")+
      sjPlot::theme_blank()+
      theme(legend.position="bottom")
## 
## -- Column specification --------------------------------------------------------
## cols(
##   Fecha = col_character(),
##   Region = col_character(),
##   `Nuevo Confirmado` = col_double(),
##   `Nuevo Muerte` = col_double(),
##   `Nuevo Recuperado` = col_double(),
##   `Acum Confirmado` = col_double(),
##   `Acum Muerte` = col_double(),
##   `Acum Recuperado` = col_double(),
##   `Casos UCI` = col_double()
## )
covid19_chile_resto <-
    readr::read_csv("https://raw.githubusercontent.com/ivanMSC/COVID19_Chile/master/covid19_chile.csv") %>% 
    dplyr::select(Fecha, Region, `Nuevo Confirmado`) %>%
    dplyr::mutate(Fecha=as.Date.character(Fecha,tryFormats = "%d-%m-%Y")) %>% 
    dplyr::filter(!Region %in% c("Valparaíso","Coquimbo")) %>% 
    dplyr::filter(!Region %in% c("Metropolitana")) %>% 
    dplyr::group_by(Fecha)%>% 
    dplyr::summarise(mdn=median(`Nuevo Confirmado`,na.rm=T),
                     q1=quantile(`Nuevo Confirmado`,.25,na.rm=T),
                     q3=quantile(`Nuevo Confirmado`,.75,na.rm=T))
## 
## -- Column specification --------------------------------------------------------
## cols(
##   Fecha = col_character(),
##   Region = col_character(),
##   `Nuevo Confirmado` = col_double(),
##   `Nuevo Muerte` = col_double(),
##   `Nuevo Recuperado` = col_double(),
##   `Acum Confirmado` = col_double(),
##   `Acum Muerte` = col_double(),
##   `Acum Recuperado` = col_double(),
##   `Casos UCI` = col_double()
## )
covid19_chile_coq_val+ 
    geom_ribbon(data= covid19_chile_resto, aes(y= mdn, x= Fecha, ymin = q1, ymax = q3), fill = "grey70", color="grey70", alpha=.4)+
    geom_line(data= covid19_chile_resto, aes(x= Fecha, y= mdn, color="black"))+
    scale_color_manual(name= "Región", values = c("black","red","blue"), 
                       label=c("Mediana (resto regiones\nexcluyendo Metropolitana)",
                               "Coquimbo",
                               "Valparaíso"))

  • ¿No se ve más ordenado?, lo valorarán cuando trabajen con condiciones anidadas (IFs concatenados)

Otra forma de escribirlo

Factorizamos lo que tienen en común las bases de datos, y para cada capa filtramos los datos de interés.

covid19_chile <-
    readr::read_csv("https://raw.githubusercontent.com/ivanMSC/COVID19_Chile/master/covid19_chile.csv") %>% 
    dplyr::select(Fecha, Region, `Nuevo Confirmado`) %>%
    dplyr::mutate(Fecha=as.Date.character(Fecha,tryFormats = "%d-%m-%Y"))
## 
## -- Column specification --------------------------------------------------------
## cols(
##   Fecha = col_character(),
##   Region = col_character(),
##   `Nuevo Confirmado` = col_double(),
##   `Nuevo Muerte` = col_double(),
##   `Nuevo Recuperado` = col_double(),
##   `Acum Confirmado` = col_double(),
##   `Acum Muerte` = col_double(),
##   `Acum Recuperado` = col_double(),
##   `Casos UCI` = col_double()
## )
ggplot()+
      geom_line(data=covid19_chile %>% 
                  dplyr::filter(Region %in% c("Coquimbo","Valparaíso")),
                aes(x=as.Date(Fecha), y=`Nuevo Confirmado`, color=Region))+
      scale_x_date()+
      xlab("Fecha")+
      sjPlot::theme_blank()+
      theme(legend.position="bottom")+
      geom_ribbon(data= covid19_chile %>% 
                  dplyr::filter(!Region %in% c("Coquimbo","Valparaíso","Metropolitana")) %>% 
                  dplyr::group_by(Fecha)%>% 
                  dplyr::summarise(mdn=median(`Nuevo Confirmado`,na.rm=T),
                     q1=quantile(`Nuevo Confirmado`,.25,na.rm=T),
                     q3=quantile(`Nuevo Confirmado`,.75,na.rm=T)), 
                aes(y= mdn, x= Fecha, ymin = q1, ymax = q3), fill = "grey70", color="grey70", alpha=.4)+
    geom_line(data= covid19_chile_resto, aes(x= Fecha, y= mdn,color="black"))+
    scale_color_manual(name= "Región", values = c("black","red","blue"), 
                       label=c("Mediana (resto regiones\nexcluyendo Metropolitana)",
                               "Coquimbo",
                               "Valparaíso"))

Facet Wrap por Región

Muy complejo y todo, pero no entendí nada los gráficos anteriores. Si queremos separar los datos por región:

ggplot()+
      geom_line(data=covid19_chile %>% 
                  dplyr::filter(Region %in% c("Coquimbo","Valparaíso")),
                aes(x=as.Date(Fecha), y=`Nuevo Confirmado`, color=Region))+
        geom_ribbon(data= covid19_chile %>% 
                  dplyr::filter(!Region %in% c("Coquimbo","Valparaíso","Metropolitana")) %>% 
                  dplyr::group_by(Fecha)%>% 
                  dplyr::summarise(mdn=median(`Nuevo Confirmado`,na.rm=T),
                     q1=quantile(`Nuevo Confirmado`,.25,na.rm=T),
                     q3=quantile(`Nuevo Confirmado`,.75,na.rm=T)), 
                aes(y= mdn, x= Fecha, ymin = q1, ymax = q3), fill = "grey70", color="grey70", alpha=.4)+
      scale_x_date()+
      xlab("Fecha")+
      sjPlot::theme_blank()+
      theme(legend.position="bottom")+
      facet_wrap(~Region, ncol = 1)

Ejercicio

  • Ocupe la misma base de datos.

  • Seleccione las muertes diarias notificadas (Nuevo Muerte)

  • Haga un gráfico con Valparaíso en comparación a la mediana y el rango intercuartil de los conteos diarios del resto de las regiones en conjunto.

  • Nótese la calidad de los datos provistos (ej., 2020-06-08, ¿el 50% obtuvo menos 1,5 muertes?)

Vacunaciones COVID-19(1)

  • Existen antecedentes de que gran parte de las parejas en Coquimbo y Valparaíso asistió a pubs y restaurantes dado que estaban en mayor proporción abiertos. El resto, recurrió a otros medios de contacto
  • Esto no ocurrió en regiones, quienes se encontraban en cuarentena estricta
  • Dado que algunos expertos señalan que la transmisibilidad se reduce con la primera dosis, mayores vacunaciones debiesen confundir las asociaciones
  • ¿El 14 de febrero en Coquimbo y Valparaíso tuvo un efecto en el aumento de casos (controlando por vacunaciones)?
  • Debemos incorporar las vacunaciones a la base de datos de nuevos casos
  • Otro paquete de la familia Tidyverse es “tidyr”
  • Similar a las funciones utilizadas por reshape (cambiar el formato largo y ancho de los datos)
  • También mostraremos como ejemplo los join.
  • Permiten combinar bases de datos
    • Left
    • Inner
    • Right
    • Otros (ej., probabilístico)

Inferencia Causal

  • “scientific generalization is a broader question than mathematical description.”
  • No importa tanto el estimador puntual, si no su varianza (intervalos de confianza, márgenes de error, poder estadístico, sesgo o error sistemático).
  • Hay muchas fuentes de sesgo que hay que atender, e incluso en distintas etapas del proceso de investigación.
  • Aproximación cuasi-experimental: Distinguir efecto causal y resultado bajo observación pasiva. No puedo aislar ni manipular variables.
  • En algunos casos se utiliza para predecir, pero es distinto poner énfasis en el predictor (B o beta, inferencia causal) que en el valor predicho (Y) (Ver diferencias entre asociación predivtiva y una relación causal).
Figura. Varianza vs. Sesgo (Fuente: Singh, 2018, Mayo 21)

Figura. Varianza vs. Sesgo (Fuente: Singh, 2018, Mayo 21)

Panel

  • Datos de panel: unidad de observación repetida en el tiempo
  • Balanceados vs. Desbalanceados
  • Las unidades tienen patrones de comportamiento individuales
library(DiagrammeR) 
# Nodes
 #node [shape = box]
 # S [label = 'Matched\n(S=1)',fontsize=7]
 # C [label = 'Not censored\n(C=0)',fontsize=7]
gr1<-
DiagrammeR::grViz("
digraph a_nice_graph  {

# Nodes
  node [shape = plaintext]
  a [label = 'Confirmados \nPre 14 Feb\nZ(0)',fontsize=10]
  b [label = 'Coquimbo y\nValparaíso\n(A)',fontsize=10]
  c [label = 'Confirmados \nPost 14 Feb\nZ(1)',fontsize=10]
  d [label = 'Confusores no-observados\n(U)',fontsize=10]
  e [label = 'Confusores medidos\nPre 14 Feb\n(C)',fontsize=10]
  blank [label = '', width = 0.01, height = 0.01]
  g [label = '', fontsize=10, width = 0.001, height = 0.001, color=White]

# Edges
  edge [color = black]

  b -> c 
  b -> blank [rankdir = TD; arrowhead = none; color = white]
  blank -> d [arrowhead = none; color = white]
  a -> c [dir= both; style = dashed]
  d -> { b a c }
  e -> { b a c }
#  d -> g [ dir = none,  color = 'white',fontcolor = white,shape=none, width=0, height=0];

  { rank = same; rankdir = LR; b; a; c }
  { rankdir = TD; a; d; e }

# Graph
  graph [overlap = true]
}")
gr1

Gráfico Acíclico Dirigido, DiD (Diferencia-en-Diferencias)

#a-A=Trat, b-Z0, c-Z1, d-U, e-C

#Fig. 2. Directed acyclic graph depicting the causal association between the treatment A, pre-exposure outcome Z(0), post-exposure outcome Z(1), measured pre-exposure confounders C, and unmeasured confounders U.

Vacunaciones COVID-19(2)

  • Se estandarizaron las variables de acuerdo a la población de la región, por cada 100,000 hab.
  • Se generó un modelo de efectos fijos para capturar eventuales cambios en coquimbo en comparación a otras regiones a partir del 14 de febrero
  • ¿Fue la forma correcta de incorporar las bases de datos? ¿Se perdió información?
library(tidyr)
library(lmtest)
library(plm)
library(broom)
library(kableExtra)

#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#Generación y compilación de las bases de datos
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:

covid19_chile_vacuna_vs_confirmados <-
    readr::read_csv("https://raw.githubusercontent.com/MinCiencia/Datos-COVID19/master/output/producto76/vacunacion.csv") %>%
  #Formato largo de fechas
    tidyr::pivot_longer(cols=starts_with("20"), names_to="fecha", values_to="recuento") %>% 
  #Formato ancho por tipo de dosis (diferenciando primera y segunda en columnas distintas)
    tidyr::pivot_wider(names_from=starts_with("Dosis"), values_from="recuento") %>% 
    dplyr::mutate(fecha=as.Date.character(fecha,tryFormats = "%Y-%m-%d")) %>% 
  # Unir con la base de datos de casos por Región
    dplyr::left_join(covid19_chile, by=c("fecha"="Fecha", "Region"="Region")) %>% 
    dplyr::filter(Region!="Total") %>% 
    dplyr::select(-Segunda) %>% 
    dplyr::filter(!is.na(Primera)) %>% 
    dplyr::mutate(met=factor(ifelse(Region %in% c("Coquimbo", "Valparaíso"),1,0))) %>% 
    dplyr::mutate(txtime=factor(ifelse(fecha>="2021-02-14",1,0))) %>% 
    dplyr::rename("conf"="Nuevo Confirmado")
## 
## -- Column specification --------------------------------------------------------
## cols(
##   .default = col_double(),
##   Region = col_character(),
##   Dosis = col_character()
## )
## i Use `spec()` for the full column specifications.
#Base de datos de PCR's y UCI's
covid19_chile_vacuna_vs_confirmados %>%     
  dplyr::left_join(readr::read_csv("https://raw.githubusercontent.com/MinCiencia/Datos-COVID19/master/output/producto7/PCR_std.csv"), 
                   #Los criterios para parear es la Región(i) y la fecha(t) 
                   by=c("Region", "fecha")) %>% 
  dplyr::left_join(readr::read_csv("https://raw.githubusercontent.com/MinCiencia/Datos-COVID19/master/output/producto8/UCI_std.csv") %>% 
                     #sacar variables redundantes: código de región, población y número
                     dplyr::select(-`Codigo region`,-Poblacion), 
                   by=c("Region", "fecha")) %>% 
  dplyr::mutate(conf_100mil=(conf/Poblacion)*100000) %>% 
  dplyr::mutate(Primera_100mil= (Primera/Poblacion)*100000) %>% 
  dplyr::mutate(PCR_100mil= (numero.x/Poblacion)*100000) %>% 
  dplyr::mutate(Hosp_UCI_100mil= (numero.y/Poblacion)*100000) %>% 
  dplyr::mutate(did=factor(as.numeric(met)*as.numeric(txtime)),
                did=factor(ifelse(did==4,1,0))) %>% 
  as.data.frame() %>% 
  #Para definir la base de datos
  assign("covid19_chile_vacuna_vs_confirmados_std_pcr_uci",., envir=.GlobalEnv)
## 
## -- Column specification --------------------------------------------------------
## cols(
##   Region = col_character(),
##   `Codigo region` = col_character(),
##   Poblacion = col_double(),
##   fecha = col_date(format = ""),
##   numero = col_double()
## )
## 
## -- Column specification --------------------------------------------------------
## cols(
##   Region = col_character(),
##   `Codigo region` = col_double(),
##   Poblacion = col_double(),
##   fecha = col_date(format = ""),
##   numero = col_double()
## )
paste0("La fecha mínima pareada fue ",min(covid19_chile_vacuna_vs_confirmados_std_pcr_uci$fecha),";\n", "La máxima fecha disponible fue ",max(covid19_chile_vacuna_vs_confirmados_std_pcr_uci$fecha))    
## [1] "La fecha mínima pareada fue 2020-12-24;\nLa máxima fecha disponible fue 2021-02-28"
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
# Modelo de efectos fijos
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:

did.reg_basico <- plm(conf_100mil ~ txtime + met + txtime:met, 
               data = covid19_chile_vacuna_vs_confirmados_std_pcr_uci,  index=c("Region", "fecha"), model = "within")
did.reg <- plm(conf_100mil ~ txtime + met + txtime:met + Primera_100mil + PCR_100mil + Hosp_UCI_100mil, 
               data = covid19_chile_vacuna_vs_confirmados_std_pcr_uci,  index=c("Region", "fecha"), model = "within")
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:


modelo_simple<-
parameters::model_parameters(did.reg_basico) %>% 
  data.frame() %>% 
  #Cambiamos los nombres de las variables
  dplyr::mutate(Parameter= dplyr::case_when(Parameter=="txtime1"~"14 de Febrero hasta el día de obtención",
                                       Parameter=="Primera_100mil"~"Primera dosis por cada 100 mil hab.",
                                       Parameter=="PCR_100mil"~"PCR's por cada 100 mil hab.",
                                       Parameter=="Hosp_UCI_100mil"~"Hospitalizaciones por cada 100 mil hab.",
                                       Parameter=="txtime1:met1"~"Coquimbo y Valparaíso x 14 de Febrero"
                                       )) %>% 
  dplyr::mutate_at(vars("Coefficient","SE","CI_low", "CI_high","t", "df_error"), funs(round(.,2))) %>% 
  #Aproximamos a 3, de manera que los interesados puedan ver el valor p con 3 decimales
  dplyr::mutate_at(vars("p"), funs(round(.,3))) %>% 
  dplyr::mutate(p=ifelse(p<0.001,"<0.001",as.character(p))) %>% 
  dplyr::mutate(CI_95=paste0(sprintf("%3.2f", CI_low),", ", sprintf("%3.2f", CI_high))) %>% 
  dplyr::select(-CI_low, -CI_high, -df_error) 

modelo_final<-
parameters::model_parameters(did.reg) %>% 
  data.frame() %>% 
  #Cambiamos los nombres de las variables
  dplyr::mutate(Parameter= dplyr::case_when(Parameter=="txtime1"~"14 de Febrero hasta el día de obtención",
                                       Parameter=="Primera_100mil"~"Primera dosis por cada 100 mil hab.",
                                       Parameter=="PCR_100mil"~"PCR's por cada 100 mil hab.",
                                       Parameter=="Hosp_UCI_100mil"~"Hospitalizaciones por cada 100 mil hab.",
                                       Parameter=="txtime1:met1"~"Coquimbo y Valparaíso x 14 de Febrero"
                                       )) %>% 
  dplyr::mutate_at(vars("Coefficient","SE","CI_low", "CI_high","t", "df_error"), funs(round(.,2))) %>% 
  #Aproximamos a 3, de manera que los interesados puedan ver el valor p con 3 decimales
  dplyr::mutate_at(vars("p"), funs(round(.,3))) %>% 
  dplyr::mutate(p=ifelse(p<0.001,"<0.001",as.character(p))) %>% 
  dplyr::mutate(CI_95=paste0(sprintf("%3.2f", CI_low),", ", sprintf("%3.2f", CI_high))) %>% 
  dplyr::select(-CI_low, -CI_high, -df_error) 

modelos_combinados<-
modelo_simple %>% 
  #mezclamos las bases de datos sin perder observaciones de ambas, y aquellas filas que calzan, quedarán unidas.
  dplyr::full_join(modelo_final,by="Parameter","_mod_final") %>% 
  dplyr::select(Parameter, Coefficient.x, CI_95.x, p.x, Coefficient.y, CI_95.y, p.y)

#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
options(knitr.kable.NA = '')

modelos_combinados %>% 
  #generar tabla en formato html, definiendo título y nombre de columnas
  knitr::kable(.,format = "html", format.args = list(decimal.mark = ".", big.mark = ","),
               caption = paste0("Tabla 1. Coeficientes Diferencia en Diferencias"),
               col.names = c("Variables Independientes","Coeficiente", "IC 95%","Sig","Coeficiente", "IC 95%","Sig"),
align =c('l',rep('c', 101)))%>%
  #Añadimos otra fila como encabezado
  kableExtra::add_header_above(c("","Modelo simple"=3, "Modelo con covariables"=3)) %>% 
  #destaca la quinta fila correspondiente a la intervención para el tratado
  kableExtra::row_spec(2,bold=T,hline_after = T) %>% 
  kableExtra::add_footnote(c("Nota. IC 95%= Intervalos de Confianza al 95%",paste0("Datos obtenidos al ", format(Sys.time(), '%d %B, %Y'))), notation = "none")%>%
  kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 13)
Tabla 1. Coeficientes Diferencia en Diferencias
Modelo simple
Modelo con covariables
Variables Independientes Coeficiente IC 95% Sig Coeficiente IC 95% Sig
14 de Febrero hasta el día de obtención -3.00 -4.62, -1.37 <0.001 -0.08 -2.56, 2.41 0.952
Coquimbo y Valparaíso x 14 de Febrero 8.09 3.49, 12.69 0.001 4.07 0.81, 7.34 0.015
Primera dosis por cada 100 mil hab. 0.00 -0.00, -0.00 0.007
PCR’s por cada 100 mil hab. 0.07 0.06, 0.07 <0.001
Hospitalizaciones por cada 100 mil hab. 0.82 0.56, 1.09 <0.001
Nota. IC 95%= Intervalos de Confianza al 95%
Datos obtenidos al 01 marzo, 2021


  • Estamos viendo si las diferencias en casos confirmados entre el periodo post-intervención con el periodo anterior a la intervención (14 feb), es distinto en estas mismas regiones tratadas, respecto las diferencias pre-post de las regiones anteriores.


se <- function(x) sqrt(var(x)/length(x))

Produc_means <- covid19_chile_vacuna_vs_confirmados_std_pcr_uci %>% 
  dplyr::group_by(Region) %>% 
  transmute(y_median = median(conf_100mil),
            y = conf_100mil, 
            fecha = fecha) %>% 
  dplyr::ungroup() %>% 
  dplyr::mutate(y_pred = predict(did.reg) + y_median) %>% 
  dplyr::select(-y_median)

#OTRA FORMA DE ESTIMAR LOS PREDICTORES
#https://stackoverflow.com/questions/28547614/prediction-with-plm-method
pmodel_fixed <-pmodel.response(did.reg,model='fixed') 


Produc_means %>% 
  gather(type, value, y, y_pred) %>% 
  dplyr::mutate(tx=factor(ifelse(Region %in% c("Coquimbo", "Valparaíso"),1,0))) %>% 
  dplyr::group_by(tx,type,fecha) %>% 
  dplyr::mutate(median=median(value), se=se(value)) %>%
  dplyr::ungroup() %>% 
  ggplot(aes(x = fecha, y = median, linetype = type))+
  geom_line() +
  geom_ribbon(aes(y= median, x= fecha, ymin = median-se, ymax = median+se), fill = "grey70", color="grey70", alpha=.4)+
  scale_linetype_manual(name="Valores",values=c("solid","dashed"), labels=c("Actual", "Estimado"))+
  facet_wrap(~tx, labeller=as_labeller(c(`0` = "Otras Regiones", `1` = "Coquimbo y\nValparaíso")))+
  geom_vline(xintercept = as.Date("2021-02-14"),linetype="longdash")+
  theme_blank()+
  ggtitle("Prediciones y Medianas por fecha de nuevos casos confirmados")+
  labs(caption="Nota. Área gris representa el error estandar por región;\nLínea vertical= 14 de Febrero",
       x="Fecha",
       y="Mediana Casos Confirmados")

library(gplots)
par(mfrow=c(1,2))
plotmeans(conf_100mil ~ Region, main = "Heterogeneidad entre Regiones", ylab="Casos Confirmados (x 100,000)", data = covid19_chile_vacuna_vs_confirmados_std_pcr_uci, n.label=F)
plotmeans(conf_100mil ~ fecha, main = "Heterogeneidad entre Fechas", ylab="", data = covid19_chile_vacuna_vs_confirmados_std_pcr_uci, n.label=F)

Vacunaciones COVID-19(3)-Limitaciones

  • A partir del gráfico no es tan claro distinguir si las diferencias se deben al 14 de febrero
  • Sí es claro que las variables de control jugaron un efecto que potenció las diferencias pre-post en las comunas tratadas, en comparación al resto de regiones.
  • ¿Cuál es la distribución de los datos?
  • Las covariables podrían verse afectadas por la intervención (14 feb), en particular en las regiones tratadas
  • Adicionalmente, la capacidad de predicción del modelo se aprecia baja
  • La diferencia en cada punto de tiempo no es constante en el periodo pre-intervención
  • Sólo dos tratados
  • ¿Y la heterogeneidad intraregional? (Comunas distintas)
  • Spillover, elementos geográficos no capturados
  • Igual no son del todo comparables las series (ausencia soporte común)
  • Hay patrones temporales que no son capturados (ciclos, tendencias, estaciones, retrasos) (pbgtest())
  • Las primeras vacunas siguen un patrón muy distinto
#Creamos una función para capturar los intervalos de confianza 
cl_quantile <- function(x, q = c(0.025, 0.975), na.rm = TRUE){
  dat <- data.frame(ymin = quantile(x, probs = q[1], na.rm = na.rm),
                    ymax = quantile(x, probs = q[2], na.rm = na.rm))
  return(dat)
}

ggplot(covid19_chile_vacuna_vs_confirmados_std_pcr_uci, aes(fecha, conf_100mil, color = met)) +
    stat_summary(geom = 'line', fun.y = median) + #como línea, la mediana
    stat_summary(geom = "ribbon", fun.data = cl_quantile, alpha = 0.3, aes(fill=met))+ #ribbon, el área es la función que definimos: 95% IC. ## todo lo que no sea fijo, debe ir en "aes" 
    geom_vline(xintercept = as.Date("2021-02-14")) + #Tiempo de intervención
    scale_color_manual(name="Grupo",labels = c("Otras\nRegiones", "Coquimbo y\n Valparaíso"), values=c("red","blue"))+ #Coquimbo y Valparaíso versus el resto
      scale_fill_manual(name="Grupo",labels = c("Otras\nRegiones", "Coquimbo y\n Valparaíso"), values=c("red","blue"))+ #Si asignan los mismos valores, se mezclarán ambas especificaciones: color y relleno del área
    sjPlot::theme_sjplot2()+
  labs(caption="Nota. Área roja: intervalos de confianza al 95%;\nLínea vertical: 14 de febrero 2021",
       y="Casos confirmados (x 100.000 hab.)",
       x="Fecha")

  • ¿Soporte común?
library(tableone)
kableone <- function(x, ...) {
  capture.output(x <- print(x,...))
  knitr::kable(x,
               format= "html", 
               caption = 'Tabla 2. Descriptivos Variables de Interés, antes de 14 de Febrero',
               col.names=c("Resto de\nlas Regiones", "Coquimbo y\nValparaíso","Sig", "Prueba","SMD"), 
               format.args= list(decimal.mark= ".", big.mark= ","))
}


#attr(BD_2020_input3$latdec_dic, "label") <- 'Low Decisional Latitude (Dichotomized)(JCQ)'
#attr(BD_2020_input3$democ_dic, "label") <- 'High Job Demands (Dichotomized)(JCQ)'

tab<-
   CreateTableOne(vars = c("conf_100mil","Primera_100mil","PCR_100mil","Hosp_UCI_100mil"), 
                     data= covid19_chile_vacuna_vs_confirmados_std_pcr_uci %>% dplyr::filter(txtime==0), 
                     strata = "met", 
                     includeNA =F,
                      smd=T)

kableone(tab, nonnormal= c("conf_100mil","Primera_100mil","PCR_100mil","Hosp_UCI_100mil"),
                       test=T, varLabels=T,noSpaces=T, printToggle=T, dropEqual=F,smd=T) %>% 
   kableExtra::kable_styling(bootstrap_options = c("striped", "hover","condensed"),font_size= 12) %>%
    kableExtra::row_spec(1, bold = T) %>%
    kableExtra::add_footnote(c("Nota. SMD= Diferencia de Medias Estandarizada"), notation = "none")%>%
    kableExtra::scroll_box(width = "100%", height = "200px")
Tabla 2. Descriptivos Variables de Interés, antes de 14 de Febrero
Resto de las Regiones Coquimbo y Valparaíso Sig Prueba SMD
n 728 104
conf_100mil (median [IQR]) 25.84 [16.50, 37.51] 11.53 [7.98, 14.04] <0.001 nonnorm 1.445
Primera_100mil (median [IQR]) 71.63 [0.00, 1211.60] 74.18 [0.00, 1048.56] 0.213 nonnorm 0.090
PCR_100mil (median [IQR]) 303.10 [216.48, 407.29] 182.75 [119.65, 230.35] <0.001 nonnorm 1.234
Hosp_UCI_100mil (median [IQR]) 5.61 [3.84, 7.30] 3.62 [2.69, 5.31] <0.001 nonnorm 0.874
Nota. SMD= Diferencia de Medias Estandarizada

Fuentes

  • Amrhein, V., Trafimow, D., & Greenland, S. (2019). Inferential Statistics as Descriptive Statistics: There Is No Replication Crisis if We Don’t Expect Replication, The American Statistician, 73:sup1, 262-270, DOI: 10.1080/00031305.2018.1543137
  • Cunningham, S. (2020). Online version of Causal Inference: The Mixtape. https://mixtape.scunning.com/index.html
  • Mummolo, J. 150C Causal Inference. Difference in Differences. Lectures [presentation]. https://scholar.princeton.edu/sites/default/files/jmummolo/files/did_jm.pdf
  • Bilach, T. (2020, Mayo 29). The difference between DID and fixed effect model. https://stats.stackexchange.com/q/454261
  • Singh, S. (2018, Mayo 21). Understanding the Bias-Variance Tradeoff. Towards data science. https://towardsdatascience.com/understanding-the-bias-variance-tradeoff-165e6942b229
  • Sofer, T., Richardson, D., Colincino, E., Schwartz, J. & Tchetgen, E. (2015). On Simple Relations Between Difference-in-differences and Negative Outcome Control of Unobserved Confounding. Harvard University Biostatistics Working Paper Series, 194,1-29.doi:awesomedoi/20150801˙TSofer˙DID
  • Jepsen, P., Vilstrup, H. and Andersen, P.K. (2015), The clinical course of cirrhosis: The importance of multistate models and competing risks analysis. Hepatology, 62: 292-302. https://doi.org/10.1002/hep.27598
  • Illari, P., & Russo, F. (2014). Causality: Philosophical theory meets scientific practice. Oxford, UK: Oxford University Press.